공간정보분석
ASS-05
‘Regression: A First Cut’
※ 본 문서는 PC 환경에 최적화되어있습니다.
1 들어가며
1.1 사용한 라이브러리
1.1.1 데이터 손질
- tidyverse: 데이터 모델링, 가공, 시각화
- readxl: 엑셀데이터 읽기
- lubridate: 날짜처리 생산성 패키지
1.1.2 시각화
- ggplot2: 그래픽 패키지
- DT: 표 시각화
- ggformula: ggplot2에 formula 인터페이스 제공
- mosaic: 모자이크 플롯
- lattice: 다차원 데이터 시각화
1.1.3 맵핑 관련
- ggmap: 온라인 소스를 정적인 지도에 올려줌
- sf: simple features, 피처스 데이터 사용
- tmap: 지도 도식
- rgdal: 지리공간데이터 처리
- raster: 래스터 포멧의 공간자료 처리
1.1.4 분석 관련
- DescTools: 탐색적 데이터 분석에 유용한 패키지
- car: 응용 회귀 동반자(R Companion to Applied Regression)
- psych: 심리학 통계 패키지
- gvlma: 가설검정
- stargazer: 회귀 리포트 작성
2 대기질-교통량 회귀분석
2.1 자료정제
구득한 자료의 종류는 크게 두 가지이며, 미세먼지 등의 관측치가 포함된 서울시 대기질 자료와 서울시 주요 길목의 교통량 데이터입니다.
분석을 진행하기 전, 필요한 부분을 취하고 버릴 부분은 도려낸 뒤, 관리가 쉽도록 하나의 데이터프레임으로 추출하였습니다.
2.1.1 공기질 데이터 불러오기
구득한 엑셀데이터를 ’readxl’라이브러리를 사용하여 불러온 후, 측정소명이 ’로’로 끝나는 행만 가져왔습니다.
AQraw <- read_xlsx("./data/일별평균대기오염도_2018.xlsx")
AQroad <- subset(AQraw, grepl(glob2rx("*로"), 측정소명))
head(AQroad)# A tibble: 6 x 8
측정일시 측정소명 `이산화질소농도(ppm)`~ `오존농도(ppm)` `이산화탄소농도(ppm)`~
<dbl> <chr> <dbl> <dbl> <dbl>
1 20180101 강남대로 0.04 0.007 0.8
2 20180101 강변북로 0.033 0.008 0.6
3 20180101 공항대로 0.045 0.007 0.6
4 20180101 도산대로 0.037 0.009 0.6
5 20180101 동작대로 0.051 0.005 0.7
6 20180101 신촌로 0.045 0.005 1
# ... with 3 more variables: `아황산가스(ppm)` <dbl>, `미세먼지(㎍/㎥)` <dbl>,
# `초미세먼지(㎍/㎥)` <dbl>
2.1.2 관측지점 자료와의 조인
구득한 관측지점을 시각화하고자 관측지점 주소 데이터를 조인하였습니다.
AQlocation <- read_xls("./data/station_list.xls", sheet = 1)[-c(1, 2), -c(1, 4, 5)]
AQlocation[5, 1] <- "동작대로"
names(AQlocation) <- c("측정소명", "주소")
head(AQlocation)# A tibble: 6 x 2
측정소명 주소
<chr> <chr>
1 강남대로 서울특별시 서초구 강남대로 201 서초구민회관 앞 중앙차로 (양재동)
2 강변북로 서울 성동구 강변북로 257 한강사업본부 옆
3 공항대로 서울 강서구 마곡동 727-1091 마곡역 중앙차로정류장 옆
4 도산대로 서울 강남구 도산대로 104 신사역2번출구 앞
5 동작대로 서울 동작구 동작대로 144 이수역 북단 버스중앙차로
6 시흥대로 서울 금천구 독산동 996-9 한양수자인아파트 앞
# A tibble: 4,362 x 9
측정일시 측정소명 `이산화질소농도(ppm)`~ `오존농도(ppm)` `이산화탄소농도(ppm)`~
<dbl> <chr> <dbl> <dbl> <dbl>
1 20180101 강남대로 0.04 0.007 0.8
2 20180101 강변북로 0.033 0.008 0.6
3 20180101 공항대로 0.045 0.007 0.6
4 20180101 도산대로 0.037 0.009 0.6
5 20180101 동작대로 0.051 0.005 0.7
6 20180101 신촌로 0.045 0.005 1
7 20180101 영등포로 0.049 0.007 0.6
8 20180101 정릉로 0.05 0.007 0.8
9 20180101 종로 NA NA NA
10 20180101 천호대로 0.045 0.007 0.6
# ... with 4,352 more rows, and 4 more variables: `아황산가스(ppm)` <dbl>,
# `미세먼지(㎍/㎥)` <dbl>, `초미세먼지(㎍/㎥)` <dbl>, 주소 <chr>
2.1.3 조인을 위한 인덱스 생성 1
같은 위치상에 존재하는 대기질 관측소와 교통량 정보 수집 장치의 정보를 구할 수 없어, 비슷한 위치상의 것들로 대체하였습니다. 주소가 달라 조인 인덱스로 사용할 수 없었기에, 날짜와 도로 명을 합친 인덱싱을 위한 하나의 새로운 열을 생성하였습니다.
2.1.4 교통량 데이터 불러오기
구득한 데이터는 월별로 되어있어, for문을 사용하여 불러온 뒤 열을 기준으로 마지막 행 뒤에 차례로 합하였습니다.
[1] "./data/01월 서울시 교통량 조사자료.xlsx"
[2] "./data/02월 서울시 교통량 조사자료.xlsx"
[3] "./data/03월 서울시 교통량 조사자료.xlsx"
for (g in 1:12) {
assign(paste0("T", g, "raw"), read_xlsx(fi[g]))
}
Traw <- bind_rows(T1raw, T2raw, T3raw, T4raw, T5raw, T6raw, T7raw, T8raw, T9raw,
T10raw, T11raw, T12raw)
Traw# A tibble: 105,304 x 30
일자 요일 지점명 지점번호 구분 방향 `0시` `1시` `2시` `3시` `4시` `5시`
<dbl> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 2.02e7 월 성산로(금~ A-01 유입 봉원고가~ 657 617 405 279 266 365
2 2.02e7 화 성산로(금~ A-01 유입 봉원고가~ 382 256 166 171 266 751
3 2.02e7 수 성산로(금~ A-01 유입 봉원고가~ 641 492 310 238 329 691
4 2.02e7 목 성산로(금~ A-01 유입 봉원고가~ 661 483 357 279 357 644
5 2.02e7 금 성산로(금~ A-01 유입 봉원고가~ 764 524 361 282 339 672
6 2.02e7 토 성산로(금~ A-01 유입 봉원고가~ 807 594 486 412 372 521
7 2.02e7 일 성산로(금~ A-01 유입 봉원고가~ 617 478 354 264 230 316
8 2.02e7 월 성산로(금~ A-01 유입 봉원고가~ 448 278 231 205 256 715
9 2.02e7 화 성산로(금~ A-01 유입 봉원고가~ 641 467 304 256 312 679
10 2.02e7 수 성산로(금~ A-01 유입 봉원고가~ 589 451 291 262 295 569
# ... with 105,294 more rows, and 18 more variables: `6시` <dbl>, `7시` <dbl>,
# `8시` <dbl>, `9시` <dbl>, `10시` <dbl>, `11시` <dbl>, `12시` <dbl>,
# `13시` <dbl>, `14시` <dbl>, `15시` <dbl>, `16시` <dbl>, `17시` <dbl>,
# `18시` <dbl>, `19시` <dbl>, `20시` <dbl>, `21시` <dbl>, `22시` <dbl>,
# `23시` <dbl>
2.1.5 두 데이터간 상응지점 지정
위에서도 언급하였듯이, 같은 위치상의 대기질 관측소와 교통량 수집기 데이터를 구득하지 못하였기에, 명칭과 주소를 대조하여 아래와 같이 11개의 서로 비슷한 지점을 찾을 수 있었습니다.
해당 작업은 if문을 사용하여 교통량 데이터의 지점에 상응하는 지점번호(코드)를 기준으로 대조를 통해 찾은 관측소들을 교통량 데이터 데이터프레임에 추가하였습니다.
for (i in 1:length(Traw$지점번호)) {
if (Traw$지점번호[i] == "A-13") {
# 1
a[i] <- "한강대로"
} else if (Traw$지점번호[i] == "A-15") {
# 2
a[i] <- "종로"
} else if (Traw$지점번호[i] == "D-09") {
# 3
a[i] <- "정릉로"
} else if (Traw$지점번호[i] == "D-10") {
# 4
a[i] <- "화랑로"
} else if (Traw$지점번호[i] == "D-21") {
# 5
a[i] <- "공항대로"
} else if (Traw$지점번호[i] == "D-24") {
# 6
a[i] <- "시흥대로"
} else if (Traw$지점번호[i] == "D-25") {
# 7
a[i] <- "영등포로"
} else if (Traw$지점번호[i] == "D-31") {
# 8
a[i] <- "동작대로"
} else if (Traw$지점번호[i] == "D-35") {
# 9
a[i] <- "강남대로"
} else if (Traw$지점번호[i] == "D-43") {
# 10
a[i] <- "도산대로"
} else if (Traw$지점번호[i] == "F-02") {
# 11
a[i] <- "강변북로"
} else {
# else
a[i] <- "NULL"
}
}
cat_traffic <- cbind(Traw, a)
# remove null value in column a
cat1_traffic <- cat_traffic[!(cat_traffic$a == "NULL"), ]2.1.6 일별 전체 합 구하기
0시 1시 2시 3시 4시 5시 6시 7시 8시 9시 10시 11시 12시 13시 14시 15시
745 514 458 258 288 194 360 549 535 783 883 1053 1218 1360 1383 1350 1411
746 326 287 189 176 230 539 1230 1976 2223 2164 1877 1936 1714 2067 1899 1911
747 565 406 300 211 223 519 1100 1986 2190 2124 2076 2101 2026 2037 1989 2068
16시 17시 18시 19시 20시 21시 22시 23시
745 1287 1192 986 881 809 777 727 575
746 1835 1815 1809 1504 1219 1092 996 761
747 1842 1742 1729 1638 1193 1091 953 756
[ reached 'max' / getOption("max.print") -- omitted 3 rows ]
교통량 데이터가 하루 단위가 아닌 시간 단위로 되어 있고, 유출입이 나뉘어 있어, 하루 단위인 대기질 데이터와 맞춰줄 필요가 있었습니다.
위에서 먼저 시간순으로 정렬하여 지점별 유출입 행끼리 마주보게 하였습니다.
다음은 반복문을 사용하여 일별 총합을 구하는 과정입니다.
dsum <- vector(mode = "numeric", length(order_traffic$a))
for (j in 1:length(order_traffic$a)) {
tmp <- vector(mode = "numeric", 24)
for (k in 7:30) {
tmp[k - 6] <- order_traffic[j, k]
}
dsum[j] <- sum(tmp)
}
iosum_traffic <- cbind(order_traffic, dsum) 일자 요일 지점명 지점번호 구분 방향 0시 1시 2시 3시
745 20180101 월 세종대로(서울역) A-13 유입 서울역->숭례문 514 458 258 288
776 20180101 월 세종대로(서울역) A-13 유출 숭례문->서울역 806 632 345 256
4시 5시 6시 7시 8시 9시 10시 11시 12시 13시 14시 15시 16시 17시 18시 19시
745 194 360 549 535 783 883 1053 1218 1360 1383 1350 1411 1287 1192 986 881
776 257 367 461 454 691 980 972 968 1179 1251 1230 1328 1388 1334 1172 1053
20시 21시 22시 23시 a dsum
745 809 777 727 575 한강대로 19831
776 1022 944 716 608 한강대로 20414
[ reached 'max' / getOption("max.print") -- omitted 4 rows ]
위처럼 지점별로 유출입 행이 마주보는 데이터 프레임을 가지고, 일일합 열을 2개씩 합하는 반복문을 사용하여 합하였습니다.
# sum dsum(io)
sumdsum <- vector(mode = "numeric", length(iosum_traffic$dsum)/2)
for (l in 1:as.numeric(length(iosum_traffic$dsum)/2)) {
if (l == 1) {
sumdsum[l] <- iosum_traffic$dsum[l] + iosum_traffic$dsum[l + 1]
} else {
sumdsum[l] <- iosum_traffic$dsum[l * 2 - 1] + iosum_traffic$dsum[l * 2]
}
}
gyotong <- iosum_traffic[c(seq(from = 1, to = length(iosum_traffic$구분), by = 2)),
]
traffic_final <- cbind(gyotong, sumdsum)[, -c(2, 4, 5, 6, 32)] 일자 지점명 0시 1시 2시 3시 4시 5시 6시 7시 8시 9시 10시 11시
745 20180101 세종대로(서울역) 514 458 258 288 194 360 549 535 783 883 1053 1218
869 20180101 종로(종로3가역) NA NA NA NA NA NA NA NA NA NA NA NA
12시 13시 14시 15시 16시 17시 18시 19시 20시 21시 22시 23시 a
745 1360 1383 1350 1411 1287 1192 986 881 809 777 727 575 한강대로
869 NA NA NA NA NA NA NA NA NA NA NA NA 종로
sumdsum
745 40245
869 NA
[ reached 'max' / getOption("max.print") -- omitted 4 rows ]
2.1.7 조인을 위한 인덱스 생성 2
대기질 데이터와 마찬가지로 날짜와 지정한 상응지점 값을 합하여 인덱싱을 위한 새로운 열을 생성하였습니다.
2.1.8 대기질-교통량 데이터 조인
생성한 인덱스 열을 바탕으로 두 데이터를 조인하였습니다. inner-join을 사용하여 두 인덱스 열에 존재하지 않는 데이터는 자연스럽게 사라지게 하였습니다. 조인 후에 필요하지 않을 것이라 판단되는 열은 제거하였습니다.
2.2 데이터 시각화(맵핑)
이렇게 뽑아낸 데이터의 공간적 분포를 한눈에 볼 수 있도록 맵핑을 하였습니다.
2.2.1 지형자료로 변환을 위한 손질
지형자료 속성 테이블로의 변환을 위한 전처리 과정입니다.
# defining specific time to extract each observation station from time series
# data.
AQinfo <- subset(data_refined, 측정일시 == 20181231)
# traffic collector locations
trafinfo <- read_xlsx("./data/서울시 교통량 수집정보.xlsx", sheet = 2)[-c(136, 137),
-c(3, 4)] %>% cbind(read_excel("./data/서울시 교통량 수집정보.xlsx", sheet = 2,
col_types = "numeric")[-c(136, 137), c(3, 4)]) %>% rename(지점명 = "지점명칭")
trafinfo1 <- inner_join(trafinfo, AQinfo, by = "지점명")[, c(1, 2, 8, 7)]2.2.2 지오코딩
R에서 지오코딩을 하기 위하여 Google의 API를 사용하였습니다. (API 키는 개인 키라 가렸습니다.)
register_google("구글API키주세요!") # 과제 게시시 해당 청크 삭제할 것!
AQinfo1 <- as.data.frame(AQinfo) # tibble will not work 측정일시 측정소명 이산화질소농도(ppm) 오존농도(ppm) 이산화탄소농도(ppm)
1 20181231 강남대로 0.063 0.004 1.2
2 20181231 강변북로 0.057 0.004 0.6
3 20181231 도산대로 0.037 0.006 1.3
4 20181231 동작대로 0.030 0.004 0.9
5 20181231 시흥대로 0.056 0.006 1.0
아황산가스(ppm) 미세먼지(㎍/㎥) 초미세먼지(㎍/㎥)
1 0.005 52 31
2 0.005 50 29
3 0.004 48 27
4 0.008 48 26
5 0.005 44 28
주소
1 서울특별시 서초구 강남대로 201 서초구민회관 앞 중앙차로 (양재동)
2 서울 성동구 강변북로 257 한강사업본부 옆
3 서울 강남구 도산대로 104 신사역2번출구 앞
4 서울 동작구 동작대로 144 이수역 북단 버스중앙차로
5 서울 금천구 독산동 996-9 한양수자인아파트 앞
지점명 sumdsum lon lat
1 강남대로(강남역-신분당) 80928 127.0360 37.48183
2 강변북로 230511 127.0416 37.53896
3 강남대로(신사역) 94804 127.0201 37.51605
4 동작대로(총신대입구역) 107203 126.9819 37.48649
5 시흥대로(시흥IC) 78455 126.8996 37.47557
[ reached 'max' / getOption("max.print") -- omitted 4 rows ]
2.2.3 좌표계 정의
좌표계를 ’WGS1984’로 정의하기 위해 TA 지오코딩 세션에서 구득한 ’tl_emd_seoul_4326.shp’의 좌표계 정보를 사용하였습니다.
# OBTAIN CRS from existing shapefile
sigungu_shp <- readOGR(dsn = "./data", layer = "tl_emd._seoul_4326", encoding = "UTF-8",
stringsAsFactors = FALSE)OGR data source with driver: ESRI Shapefile
Source: "C:\Users\dydgn\OneDrive - 공주대학교\바탕 화면\Analytical_Methods_for_Spatial_Information_Assignment\Zulu-05\data", layer: "tl_emd._seoul_4326"
with 467 features
It has 3 fields
# define projection
coor <- SpatialPointsDataFrame(gc[, c(12:13)], gc, proj4string = crs)
tcoor <- SpatialPointsDataFrame(trafinfo1[, c(3:4)], trafinfo1, proj4string = crs)CRS arguments:
+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
CRS arguments:
+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
잘 적용이 되었는지 plot을 해보았습니다.
2.2.4 Shape 저장 및 지도 도식
추후 사용 가능성을 염두에 두고 셰이프로 저장하고
# writeOGR(coor, paste0(getwd(),'/data/ex_shapefiles'), 'AQstations',
# driver='ESRI Shapefile', layer_options = 'encoding=UTF-8') writeOGR(tcoor,
# paste0(getwd(),'/data/ex_shapefiles'), 'traffic', driver='ESRI Shapefile',
# layer_options = 'encoding=UTF-8')위에서 정제한 데이터로 ’tmap’라이브러리를 사용하여 도식하였습니다.
# echo map
tmap_mode("view")
tm_shape(coor) + tm_markers("측정소명", id = "측정소명", size = 0.5, shape = tmap_icons("./fig/lung.png"),
legend.col.show = FALSE) + tm_shape(tcoor) + tm_markers("지점명", id = "지점명",
size = 0.5, shape = tmap_icons("./fig/car.png"), legend.col.show = FALSE)
2.3 단순회귀분석
2.3.1 장소 선택
서울의 여러 관측소에서 나타난 일별 미세먼지 농도, 이산화질소 농도, 교통량 데이터 중 강변북로, 종로, 동작대로 관측소 데이터의 상관관계를 회귀분석을 통해 알아보기로 했습니다.
2.3.2 산점도 행렬 분석
강변북로, 종로, 동작대로의 여러 공기의 질 관련 데이터 중 산점도 행렬 분석(pairs.panels)을 이용해 각 변수와 교통량의 관계와 회귀분석 과정에서 나타나는 잔차의 분포에 대해서 알아봤습니다.
잔차는 그래프를 보면 알 수 있듯이, 정확하지는 않지만 정규분포와 유사한 양상으로 나타났습니다.
따라서, 선형관계라는 가정이 적절하다는 것을 확인할 수 있습니다.
2.3.3 밀도함수
강변북로, 종로, 동작대로 세 장소 미세먼지 농도의 측정값을 바탕으로 밀도함수를 그려봤습니다.
2.3.4 미세먼지-교통량 회귀분석
회귀분석을 통해 강변북로, 종로, 강남대로 세 장소의 교통량과 미세먼지 농도의 상관관계에 대해 알아봤습니다.
2.3.4.1 강변북로
gf_point(gangbyun$`미세먼지(㎍/㎥)` ~ gangbyun$sumdsum, data = gangbyun) %>% gf_lm(interval = "prediction",
fill = "red") %>% gf_lm(interval = "confidence", fill = "navy") + geom_text(x = 215000,
y = 110, label = "y = 0.0008711x -167.4") + geom_text(x = 215000, y = 102, label = "p-value = 7.707e-09, R²= 0.1064") +
labs(x = "통행량(대/일)", y = "미세먼지(㎍/㎥)")2.3.4.2 종로
gf_point(jongno$`미세먼지(㎍/㎥)` ~ jongno$sumdsum, data = jongno) %>% gf_lm(interval = "prediction",
fill = "red") %>% gf_lm(interval = "confidence", fill = "navy") + geom_text(x = 35000,
y = 100, label = "y = 0.001892x -51.91") + geom_text(x = 35000, y = 92, label = "p-value = 0.0001899, R²= 0.09567") +
labs(x = "통행량(대/일)", y = "미세먼지(㎍/㎥)")2.3.4.3 동작대로
gf_point(dongjak$`미세먼지(㎍/㎥)` ~ dongjak$sumdsum, data = dongjak) %>% gf_lm(interval = "prediction",
fill = "red") %>% gf_lm(interval = "confidence", fill = "navy") + geom_text(x = 1e+05,
y = 105, label = "y = 0.001169x -85.21") + geom_text(x = 1e+05, y = 97, label = "p-value = 0.0002134, R²= 0.08916") +
labs(x = "통행량(대/일)", y = "미세먼지(㎍/㎥)")회귀분석 결과, 세 장소 모두 일차함수 식, plot에 나타나는 그래프 상으로 유의미한 양의 상관관계를 보임을 확인할 수 있습니다.
2.3.5 미세먼지-교통량 회귀분석 검증
2.3.5.1 강변북로
Call:
lm(formula = gangbyun$`미세먼지(㎍/㎥)` ~ gangbyun$sumdsum)
Residuals:
Min 1Q Median 3Q Max
-46.84 -17.64 -3.32 15.22 77.22
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.674e+02 3.557e+01 -4.707 3.85e-06 ***
gangbyun$sumdsum 8.711e-04 1.465e-04 5.946 7.71e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 23.21 on 297 degrees of freedom
Multiple R-squared: 0.1064, Adjusted R-squared: 0.1034
F-statistic: 35.35 on 1 and 297 DF, p-value: 7.707e-09
회귀분석에서의 P-value(유의확률)와 r제곱 값(결정계수)을 조사해봤습니다. P-value는 7.71e-10으로 나타났는데 이는 분석 및 관계의 유의함을 결정하는 값인 기각역 0.05보다 매우 작은 값입니다. 따라서 P값을 기준으로 볼 때, 귀무가설을 기각하고 대립가설을 채택해 유의미한 상관관계 분석으로 평가할 수 있을 것입니다.
그러나, r제곱 값(결정계수)은 0.1064로 비교적 낮게 나타났습니다.
결정계수가 낮게 나타나 문제가 조금 있을 수도 있지만, 분석 자체가 크게 문제가 되기보다는 상관관계 분석에 있어서 종속변수, 즉 설명변수가 동일한 관측 값이 많았던 것이 결정계수가 낮아지는데 일부 영향을 미친 것으로 추정됩니다.
따라서, 교통량과 미세먼지 농도의 상관관계에 대한 위의 분석은 낮은 결정계수로 인해 일부 한계가 존재하지만, 충분히 의미 있는 분석이라 평가할 수 있습니다.
2.3.5.2 종로
Call:
lm(formula = jongno$`미세먼지(㎍/㎥)` ~ jongno$sumdsum)
Residuals:
Min 1Q Median 3Q Max
-31.886 -11.990 -2.771 8.452 81.063
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.191e+01 2.227e+01 -2.330 0.02123 *
jongno$sumdsum 1.892e-03 4.934e-04 3.835 0.00019 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 18.24 on 139 degrees of freedom
Multiple R-squared: 0.09567, Adjusted R-squared: 0.08916
F-statistic: 14.7 on 1 and 139 DF, p-value: 0.0001899
회귀분석에서의 P-value(유의확률)와 r제곱 값(결정계수)을 조사해본 결과, P-value는 0.0001정도로 나타났습니다. 이 값은 분석 및 관계의 유의함을 결정하는 값인 기각역 0.05보다 매우 작게 나타난 것입니다. 따라서 P값을 기준으로 볼 때 귀무가설을 기각하고 대립가설을 채택해 유의미한 분석이라고 평가할 수 있을 것입니다.
그러나, r제곱 값(결정계수)은 0.0956로 비교적 낮게 나타났습니다.
결정계수가 낮게 나타나 문제가 조금 있을 수도 있지만, 분석 자체가 크게 문제가 되기보다는 상관관계 분석에 있어서 종속변수, 즉 설명변수가 동일한 관측 값이 많았던 것이 결정계수가 낮아지는데 일부 영향을 미친 것으로 추정됩니다.
따라서, 교통량과 미세먼지 농도의 상관관계에 대한 위의 분석은 낮은 결정계수로 인해 일부 한계가 존재하지만, 충분히 의미 있는 분석이라 평가할 수 있을 것입니다.
2.3.5.3 동작대로
Call:
lm(formula = dongjak$`미세먼지(㎍/㎥)` ~ dongjak$sumdsum)
Residuals:
Min 1Q Median 3Q Max
-45.061 -17.583 -3.873 11.453 90.802
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -8.521e+01 3.501e+01 -2.434 0.015557 *
dongjak$sumdsum 1.169e-03 3.117e-04 3.751 0.000213 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 23.86 on 284 degrees of freedom
Multiple R-squared: 0.04721, Adjusted R-squared: 0.04385
F-statistic: 14.07 on 1 and 284 DF, p-value: 0.0002134
회귀분석에서의 P-value(유의확률)와 r제곱 값(결정계수)을 조사해봤습니다. P-value는 0.0002로 나타났는데 이는 분석 및 관계의 유의함을 결정하는 값인 0.05보다 매우 작은 수치입니다. 따라서 P값을 기준으로 볼 때, 귀무가설을 기각하고 대립가설을 채택할 수 있을 것입니다.
그러나, r제곱 값(결정계수)은 0.0472로 비교적 낮게 나타났습니다.
따라서, 교통량과 미세먼지 농도의 상관관계에 대한 위의 분석은 낮은 결정계수로 인해 일부 한계가 존재하지만, 어느 정도 의미 있는 분석이라 할 수 있을 것입니다.
gf_point(gangbyun$`이산화질소농도(ppm)` ~ gangbyun$sumdsum, data = gangbyun) %>%
gf_lm(interval = "prediction", fill = "red") %>% gf_lm(interval = "confidence",
fill = "navy") + geom_text(x = 215000, y = 0.07, label = "y = 0.00000055789x -0.09588") +
geom_text(x = 215000, y = 0.065, label = "p-value = 7.178e-10, R²= 0.1172") +
labs(x = "통행량(대/일)", y = "이산화질소농도(ppm)")2.3.5.4 종로
gf_point(jongno$`이산화질소농도(ppm)` ~ jongno$sumdsum, data = jongno) %>% gf_lm(interval = "prediction",
fill = "red") %>% gf_lm(interval = "confidence", fill = "navy") + geom_text(x = 35000,
y = 0.065, label = "y = 0.000000869x -0.002781") + geom_text(x = 35000, y = 0.06,
label = "p-value = 0.01181, R²=0.03786") + labs(x = "통행량(대/일)", y = "이산화질소농도(ppm)")2.3.5.5 동작대로
gf_point(dongjak$`이산화질소농도(ppm)` ~ dongjak$sumdsum, data = dongjak) %>% gf_lm(interval = "prediction",
fill = "red") %>% gf_lm(interval = "confidence", fill = "navy") + geom_text(x = 1e+05,
y = 0.08, label = "y = 0.0000008714x -0.06077") + geom_text(x = 1e+05, y = 0.076,
label = "p-value = 1.463e-05, R²= 0.06082") + labs(x = "통행량(대/일)", y = "이산화질소농도(ppm)")마찬가지로 분석 결과, 세 장소 모두에서 일차함수 그래프 상으로 유의미한 상관관계가 나타나는 것을 관찰할 수 있을 것입니다.
2.3.6 이산화질소-교통량 회귀분석 검증
2.3.6.1 강변북로
Call:
lm(formula = gangbyun$`이산화질소농도(ppm)` ~ gangbyun$sumdsum)
Residuals:
Min 1Q Median 3Q Max
-0.037627 -0.008847 0.001063 0.009473 0.038892
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -9.588e-02 2.126e-02 -4.51 9.33e-06 ***
gangbyun$sumdsum 5.578e-07 8.756e-08 6.37 7.18e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.01387 on 297 degrees of freedom
Multiple R-squared: 0.1202, Adjusted R-squared: 0.1172
F-statistic: 40.58 on 1 and 297 DF, p-value: 7.178e-10
회귀분석에서의 P-value(유의확률)와 r제곱 값(결정계수)을 조사해봤습니다. P-value는 7.18e-10으로 나타났는데 이는 분석 및 관계의 유의함을 결정하는 값인 0.05보다 매우 작게 나타난 것입니다. P값을 기준으로 볼 때, 귀무가설은 p값이 낮아 기각하고 대립가설을 채택할 수 있을 것입니다.
그러나, r제곱 값(결정계수)은 0.1202로 비교적 낮게 나타났습니다.
결정계수가 낮게 나타나 문제가 조금 있을 수도 있지만, 분석 자체가 크게 문제가 되기보다는 상관관계 분석에 있어서 종속변수, 즉 설명변수가 동일한 관측 값이 많았던 것이 결정계수가 낮아지는데 일부 영향을 미친 것으로 추정됩니다.
따라서, 교통량과 이산화질소 농도의 상관관계에 대한 위의 분석은 낮은 결정계수로 인해 일부 한계가 존재하지만, 충분히 의미 있는 분석이라 평가할 수 있을 것입니다.
2.3.6.2 종로
Call:
lm(formula = jongno$`이산화질소농도(ppm)` ~ jongno$sumdsum)
Residuals:
Min 1Q Median 3Q Max
-0.028371 -0.009724 -0.000216 0.008137 0.032833
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.781e-03 1.538e-02 -0.181 0.8568
jongno$sumdsum 8.690e-07 3.406e-07 2.551 0.0118 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.01259 on 139 degrees of freedom
Multiple R-squared: 0.04473, Adjusted R-squared: 0.03786
F-statistic: 6.509 on 1 and 139 DF, p-value: 0.01181
회귀분석에서의 P-value(유의확률)와 r제곱 값(결정계수)을 조사해봤습니다. P-value는 0.01180으로 나타났는데 이는 분석 및 관계의 유의함을 결정하는 값인 0.05보다 작은 값을 보임을 알 수 있습니다. 따라서 P값을 기준으로 볼 때, 유의미한 분석이라고 평가할 수 있을 것입니다.
그러나, r제곱 값(결정계수)은 0.0447로 낮게 나타났습니다.
강변북로에 대해 행한 분석과 달리, 결정계수가 지나치게 낮게 나타나 종로에서는 교통량과 이산화질소 농도의 상관관계의 설명력에 대해 일부 의문이 생겼습니다.
2.3.6.3 동작대로
Call:
lm(formula = dongjak$`이산화질소농도(ppm)` ~ dongjak$sumdsum)
Residuals:
Min 1Q Median 3Q Max
-0.02460 -0.01224 -0.00399 0.01151 0.03884
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.077e-02 2.219e-02 -2.739 0.00656 **
dongjak$sumdsum 8.714e-07 1.976e-07 4.411 1.46e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.01512 on 284 degrees of freedom
Multiple R-squared: 0.06412, Adjusted R-squared: 0.06082
F-statistic: 19.46 on 1 and 284 DF, p-value: 1.463e-05
회귀분석에서의 P-value(유의확률)와 r제곱 값(결정계수)을 조사해본 결과, P-value는 1.46e-05로 나타났습니다. 이는 분석 및 관계의 유의함을 결정하는 값인 0.05보다 매우 작은 것으로, P값을 기준으로 볼 때 귀무가설을 기각하고 대립가설을 채택할 수 있을 것입니다.
그러나, r제곱 값(결정계수)은 0.0641로 비교적 낮게 나타났습니다.
따라서, 교통량과 이산화질소 농도의 상관관계에 대한 위의 분석은 낮은 결정계수로 인해 일부 한계가 존재하지만, 어느 정도 의미 있는 분석이라 평가할 수 있을 것입니다.
3 미세먼지 - 시청률 회귀분석
3.1 데이터 불러오기
미세먼지 농도 데이터는 위의 분석에서와 마찬가지로 2018년 일별 미세먼지자료를 구득하였으며,
시청률 데이터는 AGB닐슨의 2018년 일요일에 방영한 지상파 3사 예능 프로그램(SBS 런닝맨, KBS 슈퍼맨이 돌아왔다, MBC 복면가왕)의 수도권 시청률 통계자료를 구득하였습니다.
측정일시 측정소명 이산화질소농도(ppm) 오존농도(ppm)
Min. :20180101 Length:12447 Min. :0.00000 Min. :0.00000
1st Qu.:20180410 Class :character 1st Qu.:0.02000 1st Qu.:0.01000
Median :20180804 Mode :character Median :0.03100 Median :0.01800
Mean :20180713 Mean :0.03256 Mean :0.01933
3rd Qu.:20181023 3rd Qu.:0.04200 3rd Qu.:0.02600
Max. :20181231 Max. :0.13000 Max. :0.24700
NA's :232 NA's :240
이산화탄소농도(ppm) 아황산가스(ppm) 미세먼지(㎍/㎥) 초미세먼지(㎍/㎥)
Min. :0.1000 Min. :0.00100 Min. : 3.00 Min. : 1.00
1st Qu.:0.4000 1st Qu.:0.00300 1st Qu.: 24.00 1st Qu.: 12.00
Median :0.5000 Median :0.00400 Median : 36.00 Median : 20.00
Mean :0.5527 Mean :0.00448 Mean : 40.53 Mean : 22.78
3rd Qu.:0.7000 3rd Qu.:0.00500 3rd Qu.: 52.00 3rd Qu.: 30.00
Max. :2.6000 Max. :0.14800 Max. :199.00 Max. :127.00
NA's :296 NA's :247 NA's :344 NA's :338
3.2 데이터 전처리
3.2.1 미세먼지
우선, 2018년 데이터 중 일요일의 대기질 데이터를 추출하였습니다.
# Extract Sunday air quality data
start_date <- as.Date("2018-01-1")
end_date <- as.Date("2018-12-31")
df <- data.frame(seq(as.Date(start_date), as.Date(end_date), by = 1))
names(df) <- "date_18"
df$weekday <- weekdays(as.Date(df$date_18))
str(df)'data.frame': 365 obs. of 2 variables:
$ date_18: Date, format: "2018-01-01" "2018-01-02" ...
$ weekday: chr "Monday" "Tuesday" "Wednesday" "Thursday" ...
sunday <- filter(df, weekday == "Sunday")
datelist <- air$측정일시
date <- ymd(datelist)
air <- air %>% mutate(date)
sunday_air <- air %>% filter(date %in% sunday$date_18)
sunday_air <- rename(sunday_air, 측정일시 = "측정일시", 미세먼지 = "미세먼지(㎍/㎥)",
초미세먼지 = "초미세먼지(㎍/㎥)", 이산화질소 = "이산화질소농도(ppm)", 오존 = "오존농도(ppm)",
이산화탄소 = "이산화탄소농도(ppm)", 아황산가스 = "아황산가스(ppm)")
head(sunday_air)# A tibble: 6 x 9
측정일시 측정소명 이산화질소 오존 이산화탄소 아황산가스 미세먼지 초미세먼지
<dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 20180107 강남구 0.048 0.003 0.7 0.006 46 42
2 20180107 강남대로 0.065 0.003 1.1 0.006 69 31
3 20180107 강동구 0.049 0.004 0.9 0.005 57 45
4 20180107 강변북로 0.041 0.004 0.9 0.005 54 25
5 20180107 강북구 0.047 0.007 0.9 0.004 49 31
6 20180107 강서구 0.056 0.004 1 0.004 67 45
# ... with 1 more variable: date <date>
일요일 대기질 데이터를 측정일시를 기준으로 그룹화하여 결측치를 제거한 후, 미세먼지 농도의 일별 평균을 구했습니다.
# mean of fine dust
n <- sunday_air %>% group_by(측정일시) %>% summarise(n = n())
mean <- sunday_air %>% group_by(측정일시) %>% summarise(미세먼지평균 = mean(미세먼지,
na.rm = TRUE))
mean_table <- inner_join(mean, n, by = "측정일시")
summary(mean_table$미세먼지평균) Min. 1st Qu. Median Mean 3rd Qu. Max.
12.33 23.02 34.00 37.73 52.00 77.36
초미세먼지 농도의 경우도 마찬가지로 평균을 구한 뒤 합쳤습니다.
## mean of ultra fine dust
mean2 <- sunday_air %>% group_by(측정일시) %>% summarise(초미세먼지평균 = mean(초미세먼지,
na.rm = TRUE))
mean_table <- inner_join(mean2, mean_table, by = "측정일시")
summary(mean_table$초미세먼지평균) Min. 1st Qu. Median Mean 3rd Qu. Max.
5.872 11.128 20.308 21.647 29.000 50.538
# A tibble: 6 x 4
측정일시 초미세먼지평균 미세먼지평균 n
<dbl> <dbl> <dbl> <int>
1 20180107 37.9 57.5 39
2 20180114 50.5 71.1 39
3 20180121 34.0 65.9 39
4 20180128 19.2 34 39
5 20180204 21 33.6 39
6 20180211 14.7 53.3 39
측정일시 초미세먼지평균 미세먼지평균 n
Min. :20180107 Min. : 5.872 Min. :12.33 Min. :39.00
1st Qu.:20180415 1st Qu.:11.128 1st Qu.:23.02 1st Qu.:39.00
Median :20180729 Median :20.308 Median :34.00 Median :39.00
Mean :20180703 Mean :21.647 Mean :37.73 Mean :40.31
3rd Qu.:20181014 3rd Qu.:29.000 3rd Qu.:52.00 3rd Qu.:39.00
Max. :20181230 Max. :50.538 Max. :77.36 Max. :46.00
3.2.2 시청률
시청률 데이터 역시 결측치를 제거한 후 조인합니다.
# viewrate data
viewrate <- read_csv("./data/viewratedata.csv")
viewrate <- viewrate %>% mutate(total = as.numeric(viewrate$Total)) %>% dplyr::select(Date,
total)
# eliminate missing values
totalrate <- rename(viewrate, 측정일시 = "Date", total_viewrate = "total")
totalrate <- na.omit(totalrate)
head(totalrate)# A tibble: 6 x 2
측정일시 total_viewrate
<dbl> <dbl>
1 20180107 36.3
2 20180114 35.5
3 20180121 31.4
4 20180128 35.2
5 20180204 32.8
6 20180218 33.1
측정일시 total_viewrate
Min. :20180107 Min. :21.1
1st Qu.:20180418 1st Qu.:26.0
Median :20180708 Median :28.0
Mean :20199832 Mean :28.3
3rd Qu.:20181010 3rd Qu.:30.8
Max. :21080325 Max. :36.3
3.2.3 완성된 데이터
준비 완료!
측정일시 초미세먼지평균 미세먼지평균 n
Min. :20180107 Min. : 5.872 Min. :12.33 Min. :39.00
1st Qu.:20180420 1st Qu.:11.215 1st Qu.:23.24 1st Qu.:39.00
Median :20180726 Median :21.346 Median :35.01 Median :39.00
Mean :20180713 Mean :22.395 Mean :38.46 Mean :40.48
3rd Qu.:20181016 3rd Qu.:29.530 3rd Qu.:52.17 3rd Qu.:39.00
Max. :20181230 Max. :50.538 Max. :77.36 Max. :46.00
total_viewrate
Min. :21.10
1st Qu.:26.27
Median :28.35
Mean :28.61
3rd Qu.:31.23
Max. :36.30
------------------------------------------------------------------------------
Describe mean_table (tbl_df, tbl, data.frame):
data frame: 40 obs. of 5 variables
40 complete cases (100.0%)
Nr ColName Class NAs Levels
1 측정일시 numeric .
2 초미세먼지평균 numeric .
3 미세먼지평균 numeric .
4 n integer .
5 total_viewrate numeric .
------------------------------------------------------------------------------
1 - 측정일시 (numeric)
length n NAs unique 0s mean meanCI
40 40 0 = n 0 2.02e+07 2.02e+07
100.0% 0.0% 0.0% 2.02e+07
.05 .10 .25 median .75 .90 .95
2.02e+07 2.02e+07 2.02e+07 2.02e+07 2.02e+07 2.02e+07 2.02e+07
range sd vcoef mad IQR skew kurt
1'123.00 366.00 0.00 444.78 595.50 -0.19 -1.34
lowest : 2.02e+07, 2.02e+07, 2.02e+07, 2.02e+07, 2.02e+07
highest: 2.02e+07, 2.02e+07, 2.02e+07, 2.02e+07, 2.02e+07
------------------------------------------------------------------------------
2 - 초미세먼지평균 (numeric)
length n NAs unique 0s mean meanCI
40 40 0 = n 0 22.395112 18.566333
100.0% 0.0% 0.0% 26.223891
.05 .10 .25 median .75 .90 .95
6.787821 9.615824 11.215476 21.346154 29.529762 39.550343 44.500328
range sd vcoef mad IQR skew kurt
44.666667 11.971837 0.534574 14.660090 18.314286 0.591248 -0.576656
lowest : 5.871795, 6.166667, 6.820513, 7.542857, 9.846154
highest: 39.456522, 40.394737, 44.325581, 47.820513, 50.538462
------------------------------------------------------------------------------
3 - 미세먼지평균 (numeric)
length n NAs unique 0s mean meanCI
40 40 0 = n 0 38.45879 32.55534
100.0% 0.0% 0.0% 44.36224
.05 .10 .25 median .75 .90 .95
13.06855 15.46923 23.24240 35.01429 52.17391 66.16410 69.10513
range sd vcoef mad IQR skew kurt
65.02564 18.45892 0.47997 18.91139 28.93152 0.41535 -0.98042
lowest : 12.33333, 12.69444, 13.08824, 14.38462, 15.58974
highest: 65.92308, 68.33333, 69.0, 71.10256, 77.35897
------------------------------------------------------------------------------
4 - n (integer)
length n NAs unique 0s mean meanCI
40 40 0 3 0 40.48 39.59
100.0% 0.0% 0.0% 41.36
.05 .10 .25 median .75 .90 .95
39.00 39.00 39.00 39.00 39.00 46.00 46.00
range sd vcoef mad IQR skew kurt
7.00 2.78 0.07 0.00 0.00 1.29 -0.30
level freq perc cumfreq cumperc
1 39 31 77.5% 31 77.5%
2 45 4 10.0% 35 87.5%
3 46 5 12.5% 40 100.0%
heap(?): remarkable frequency (77.5%) for the mode(s) (= 39)
------------------------------------------------------------------------------
5 - total_viewrate (numeric)
length n NAs unique 0s mean meanCI
40 40 0 39 0 28.613 27.458
100.0% 0.0% 0.0% 29.767
.05 .10 .25 median .75 .90 .95
23.680 24.810 26.275 28.350 31.225 33.200 35.215
range sd vcoef mad IQR skew kurt
15.200 3.610 0.126 3.781 4.950 0.153 -0.447
lowest : 21.1, 21.4, 23.8, 24.0, 24.9
highest: 33.1, 34.1, 35.2, 35.5, 36.3
3.3 분석
3.3.1 변수 검정
분석에 앞서 변수들의 정규성을 검정하였고, 검정 결과 미세먼지와 초미세먼지의 p-value가 유의 수준을 넘지 않았으므로 귀무가설 기각에 실패하였습니다.
Shapiro-Wilk normality test
data: mean_table$미세먼지평균
W = 0.94637, p-value = 0.05691
Shapiro-Wilk normality test
data: mean_table$초미세먼지평균
W = 0.94306, p-value = 0.04389
Shapiro-Wilk normality test
data: mean_table$total_viewrate
W = 0.98326, p-value = 0.8075
따라서 대기질 데이터에 log를 취하여 정규성을 갖도록 하였습니다.
# normalize variables
mean_table <- transform(mean_table, fine_log = log(미세먼지평균), ultrafine_log = log(초미세먼지평균))그 결과는 다음과 같습니다.
Shapiro-Wilk normality test
data: mean_table$fine_log
W = 0.95723, p-value = 0.1346
Shapiro-Wilk normality test
data: mean_table$ultrafine_log
W = 0.96345, p-value = 0.2193
3.3.2 가설설정
설정한 가설은
H0 : 미세먼지 농도는 시청률에 영향을 주지 않는다.
H1 : 미세먼지 농도는 시청률에 영향을 준다.
3.3.3 상관분석
회귀식을 도출하기에 앞서, cortest를 진행하였습니다. 그 결과 상관계수는 0.3541이고, p-value는 0.02494이기 때문에 두 변수간 유의한 상관관계가 있음을 알 수 있습니다.
Pearson's product-moment correlation
data: x and y
t = 2.3348, df = 38, p-value = 0.02494
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.04797772 0.59955155
sample estimates:
cor
0.3541935
3.3.4 단일 회귀 분석
미세먼지와 시청률간의 회귀분석을 진행하였습니다.
# regression of fine dust
model1 <- lm(total_viewrate ~ fine_log, data = mean_table)
summary(model1)
Call:
lm(formula = total_viewrate ~ fine_log, data = mean_table)
Residuals:
Min 1Q Median 3Q Max
-6.8905 -1.7322 -0.4328 1.6073 6.5856
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 20.000 3.728 5.364 4.24e-06 ***
fine_log 2.443 1.046 2.335 0.0249 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.42 on 38 degrees of freedom
Multiple R-squared: 0.1255, Adjusted R-squared: 0.1024
F-statistic: 5.451 on 1 and 38 DF, p-value: 0.02494
회귀 분석을 진행한 결과,
1. F-Statistic의 p-value 값이 유의 수준인 0.05보다 작게 도출되었습니다.
- Coefficients p-value값이 0.024로 이 역시 유의 수준 0.05보다 작게 도출되었습니다.
p-value에 따르면 미세먼지 농도는 시청률에 유의한 영향을 주는 것으로 해석할 수 있습니다.
- 그러나 결정계수의 값이 0.1수준으로 낮게 나왔기 때문에 이 회귀 모형은 통계적으로 완벽히 유의하다고 할 수 없습니다.
| Dependent variable: | |
| total_viewrate | |
| fine_log | 2.443** |
| (1.046) | |
| Constant | 20.000*** |
| (3.728) | |
| Observations | 40 |
| R2 | 0.125 |
| Adjusted R2 | 0.102 |
| Residual Std. Error | 3.420 (df = 38) |
| F Statistic | 5.451** (df = 1; 38) |
| Note: | p<0.1; p<0.05; p<0.01 |
gf_point(total_viewrate ~ fine_log, data = mean_table) %>% gf_lm(interval = "confidence") +
xlab("fine dust(㎍/㎥)") + ylab("total viewrate(%)")혹시나 하는 마음에 초미세먼지 데이터로도 회귀분석을 진행하였습니다.
# regression of ultrafinedust
model2 <- lm(total_viewrate ~ ultrafine_log, data = mean_table)
summary(model2)
Call:
lm(formula = total_viewrate ~ ultrafine_log, data = mean_table)
Residuals:
Min 1Q Median 3Q Max
-5.766 -1.939 -0.319 1.860 6.590
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 20.3200 2.6966 7.535 4.69e-09 ***
ultrafine_log 2.8053 0.8954 3.133 0.00333 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.26 on 38 degrees of freedom
Multiple R-squared: 0.2053, Adjusted R-squared: 0.1844
F-statistic: 9.815 on 1 and 38 DF, p-value: 0.003327
역시나 비슷한 결과가 도출되었습니다.
gf_point(total_viewrate ~ ultrafine_log, data = mean_table) %>% gf_lm(interval = "confidence") +
xlab("ultrafine dust(㎍/㎥)") + ylab("total viewrate(%)")3.3.5 원인 찾기
plot을 통해 회귀 분석 결과를 자세히 보겠습니다.
plot 상으로는 큰 문제가 없어 보입니다. 따라서 보다 정밀하게 알아보고자
gvlma패키지를 통해 다시 살펴보면, 역시 Heteroscedasticity 검정에서 통과하지 못했음을 확인할 수 있습니다.
Call:
lm(formula = total_viewrate ~ fine_log, data = mean_table)
Residuals:
Min 1Q Median 3Q Max
-6.8905 -1.7322 -0.4328 1.6073 6.5856
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 20.000 3.728 5.364 4.24e-06 ***
fine_log 2.443 1.046 2.335 0.0249 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.42 on 38 degrees of freedom
Multiple R-squared: 0.1255, Adjusted R-squared: 0.1024
F-statistic: 5.451 on 1 and 38 DF, p-value: 0.02494
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance = 0.05
Call:
gvlma(x = model1)
Value p-value Decision
Global Stat 7.763782 0.10062 Assumptions acceptable.
Skewness 0.002025 0.96410 Assumptions acceptable.
Kurtosis 0.262686 0.60828 Assumptions acceptable.
Link Function 0.115021 0.73450 Assumptions acceptable.
Heteroscedasticity 7.384050 0.00658 Assumptions NOT satisfied!
plot을 통해 확인한 outlier를 제거한 후 다시 회귀분석을 해보았습니다.
# remove outlier
fit_table <- mean_table[-c(1, 4, 10, 16, 17, 27), ]
re <- lm(total_viewrate ~ log(미세먼지평균), data = fit_table)
# again regression
plot(re)
Call:
lm(formula = total_viewrate ~ log(미세먼지평균), data = fit_table)
Residuals:
Min 1Q Median 3Q Max
-5.2213 -1.5737 -0.2117 1.1307 4.1737
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.783 3.055 5.166 1.23e-05 ***
log(미세먼지평균) 3.653 0.865 4.223 0.000186 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.49 on 32 degrees of freedom
Multiple R-squared: 0.3579, Adjusted R-squared: 0.3378
F-statistic: 17.84 on 1 and 32 DF, p-value: 0.0001863
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance = 0.05
Call:
gvlma(x = re)
Value p-value Decision
Global Stat 3.78793 0.43547 Assumptions acceptable.
Skewness 0.04969 0.82360 Assumptions acceptable.
Kurtosis 0.33388 0.56338 Assumptions acceptable.
Link Function 0.01352 0.90745 Assumptions acceptable.
Heteroscedasticity 3.39084 0.06556 Assumptions acceptable.
Call:
lm(formula = total_viewrate ~ log(미세먼지평균), data = fit_table)
Residuals:
Min 1Q Median 3Q Max
-5.2213 -1.5737 -0.2117 1.1307 4.1737
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 15.783 3.055 5.166 1.23e-05 ***
log(미세먼지평균) 3.653 0.865 4.223 0.000186 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.49 on 32 degrees of freedom
Multiple R-squared: 0.3579, Adjusted R-squared: 0.3378
F-statistic: 17.84 on 1 and 32 DF, p-value: 0.0001863
그 결과 p-value값은 더욱 작아졌고, 위에서 문제가 되었던 R Squared 값이 0.3378정도로 증가하였음을 확인 할 수 있었습니다. 설명변수로 반응변수를 약 33% 설명할 수 있습니다.
3.4 결론
미세먼지와 시청률간의 단순 회귀 분석 결과, p-value에 따르면 미세먼지 농도는 주말 예능 시청률에 영향을 주는 것으로 생각할 수 있었습니다, 그러나 결정계수의 값이 낮다는 점에서 한계가 존재하였습니다.
또한 잔차가 등분산이 아닐 경우에는 OLS가 아닌 WLS 혹은 GLS를 사용해야하지만 이상치를 제거하는 방식으로 OLS를 재진행한점에서도 한계점을 찾아볼 수 있습니다.
4 강수량 - 미세먼지 회귀분석
4.1 데이터 불러오기
강수량 데이터는 기상청의 일별 평균 강수량 데이터를 구득하였습니다.
# rain data
dust <- read_xlsx("./data/일별평균대기오염도_2018.xlsx", sheet = 1)
rain <- read_csv("./data/percipitation.csv")
dust <- dust[, c(1, 2, 7, 8)]
DescTools::Desc(rain)------------------------------------------------------------------------------
Describe rain (spec_tbl_df, tbl_df, tbl, data.frame):
data frame: 128 obs. of 4 variables
128 complete cases (100.0%)
Nr ColName Class NAs Levels
1 지점 numeric .
2 지점명 character .
3 일시 Date .
4 일강수량(mm) numeric .
------------------------------------------------------------------------------
1 - 지점 (numeric)
length n NAs unique 0s mean meanCI
128 128 0 1 0 108.00 108.00
100.0% 0.0% 0.0% 108.00
.05 .10 .25 median .75 .90 .95
108.00 108.00 108.00 108.00 108.00 108.00 108.00
range sd vcoef mad IQR skew kurt
0.00 0.00 0.00 0.00 0.00 <NA> <NA>
level freq perc cumfreq cumperc
1 108 128 100.0% 128 100.0%
heap(?): remarkable frequency (100.0%) for the mode(s) (= 108)
------------------------------------------------------------------------------
2 - 지점명 (character - dichotomous)
length n NAs unique
128 128 0 1
100.0% 0.0%
freq perc lci.95 uci.95'
서울 128 100.0% 97.1% 100.0%
' 95%-CI Wilson
------------------------------------------------------------------------------
3 - 일시 (Date)
length n NAs unique
128 128 0 128
100.0% 0.0%
lowest : 2018-01-08, 2018-01-09, 2018-01-10, 2018-01-12
highest: 2018-12-06, 2018-12-13, 2018-12-16, 2018-12-17
Weekday:
Pearson's Chi-squared test (1-dim uniform):
X-squared = 2.8125, df = 6, p-value = 0.832
level freq perc cumfreq cumperc
1 Monday 20 15.6% 20 15.6%
2 Tuesday 21 16.4% 41 32.0%
3 Wednesday 20 15.6% 61 47.7%
4 Thursday 21 16.4% 82 64.1%
5 Friday 14 10.9% 96 75.0%
6 Saturday 15 11.7% 111 86.7%
7 Sunday 17 13.3% 128 100.0%
Months:
Pearson's Chi-squared test (1-dim uniform):
X-squared = 8.5, df = 11, p-value = 0.6679
level freq perc cumfreq cumperc
1 January 11 8.6% 11 8.6%
2 February 6 4.7% 17 13.3%
3 March 10 7.8% 27 21.1%
4 April 11 8.6% 38 29.7%
5 May 14 10.9% 52 40.6%
6 June 12 9.4% 64 50.0%
7 July 10 7.8% 74 57.8%
8 August 16 12.5% 90 70.3%
9 September 12 9.4% 102 79.7%
10 October 11 8.6% 113 88.3%
11 November 9 7.0% 122 95.3%
12 December 6 4.7% 128 100.0%
By weeks :
level freq perc cumfreq cumperc
1 2018-01-08 5 3.9% 5 3.9%
2 2018-01-15 3 2.3% 8 6.2%
3 2018-01-22 1 0.8% 9 7.0%
4 2018-01-29 3 2.3% 12 9.4%
5 2018-02-05 1 0.8% 13 10.2%
6 2018-02-12 1 0.8% 14 10.9%
7 2018-02-19 2 1.6% 16 12.5%
8 2018-02-26 4 3.1% 20 15.6%
9 2018-03-05 3 2.3% 23 18.0%
10 2018-03-12 2 1.6% 25 19.5%
11 2018-03-19 2 1.6% 27 21.1%
12 2018-03-26 0 0.0% 27 21.1%
13 2018-04-02 6 4.7% 33 25.8%
14 2018-04-09 2 1.6% 35 27.3%
15 2018-04-16 1 0.8% 36 28.1%
[ reached 'max' / getOption("max.print") -- omitted 35 rows ]
------------------------------------------------------------------------------
4 - 일강수량(mm) (numeric)
length n NAs unique 0s mean meanCI
128 128 0 54 24 10.032 6.763
100.0% 0.0% 18.8% 13.301
.05 .10 .25 median .75 .90 .95
0.000 0.000 0.200 2.000 9.750 29.900 57.775
range sd vcoef mad IQR skew kurt
96.500 18.688 1.863 2.965 9.550 2.674 7.040
lowest : 0.0 (24), 0.1 (4), 0.2 (5), 0.3 (3), 0.4 (2)
highest: 64.0, 71.5, 83.0, 83.5, 96.5
heap(?): remarkable frequency (18.8%) for the mode(s) (= 0)
## 데이터 전처리
측정지점이 종로구인 데이터만 추출하였습니다.
미세먼지 데이터와 강수량 데이터를 측정일시를 기준으로 병합하였습니다.
# join data
datelist <- jongro_dust$측정일시
jongro_dust$측정일시 <- ymd(datelist)
colnames(jongro_dust)[3] <- "미세먼지"
colnames(rain)[3] <- "측정일시"
colnames(rain)[4] <- "일강수량"
rain <- merge(rain, jongro_dust, by = "측정일시")[, c(4, 6)]
rain 일강수량 미세먼지
1 0.9 49
2 0.5 37
3 0.3 27
4 0.0 21
5 0.4 48
6 0.2 63
7 0.0 97
8 0.0 80
9 3.3 38
10 2.9 36
11 0.0 35
12 0.5 31
13 0.0 44
14 0.0 41
15 0.4 42
16 3.7 79
17 25.0 67
18 0.5 35
19 0.0 52
20 11.5 28
21 3.5 8
22 0.5 25
23 4.0 29
24 27.0 18
25 0.5 52
26 1.0 30
27 5.0 37
28 9.0 18
29 15.5 33
30 59.0 9
31 3.0 7
32 0.0 52
33 12.0 41
34 1.0 24
35 22.0 13
36 32.0 27
37 0.5 18
[ reached 'max' / getOption("max.print") -- omitted 68 rows ]
4.2 분석
4.2.1 가설 설정
설정한 가설은 다음과 같습니다
HO : 강수량은 미세먼지 농도에 영향을 주지 않는다
H1 : 강수량은 미세먼지 농도에 영향을 준다
4.2.2 회귀분석
왼쪽으로 치우쳐진 분포를 보이는 일강수량과 미세먼지를 log 변환하였습니다.
# normalize variables
rain$lograin = with(rain, log(일강수량))
rain$logdust = with(rain, log(미세먼지))
par(mfrow = c(1, 2))
densityplot(~logdust, data = rain)이후 회귀 분석을 시행하였고, 회귀 분석의 결과는 다음과 같습니다.
rain$lograin[is.infinite(rain$lograin)] <- NA
lm_rain = lm(logdust ~ lograin, data = rain)
summary(lm_rain)
Call:
lm(formula = logdust ~ lograin, data = rain)
Residuals:
Min 1Q Median 3Q Max
-1.23179 -0.44525 0.01916 0.39422 1.31442
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.32133 0.07663 43.340 < 2e-16 ***
lograin -0.13392 0.03451 -3.881 0.00021 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5811 on 82 degrees of freedom
(21 observations deleted due to missingness)
Multiple R-squared: 0.1552, Adjusted R-squared: 0.1449
F-statistic: 15.06 on 1 and 82 DF, p-value: 0.0002096
p-value 값이 0.05보다 작은 0.0002 정도이므로 회귀분석 결과 귀무가설을 기각 합니다. 즉, 일 강수량은 미세먼지 농도에 영향을 준다고 볼 수 있습니다. 그러나, 결정계수 R-squared값이 약 0.14 정도로 다소 낮다는 한계를 가집니다.